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Abstract 

Stack  and  air  system  are  the  two  most  important  components  in  the  fuel  cell  system  (FCS).  It  is  meaningful  to  study  their  properties  and 
the  trade-off  between  them.  In  this  paper,  a  modified  one-dimensional  steady- state  analytical  fuel  cell  model  is  used.  The  logarithmic  mean 
of  the  inlet  and  the  outlet  oxygen  partial  pressure  is  adopted  to  avoid  underestimating  the  effect  of  air  stoichiometry.  And  the  pressure  drop 
model  in  the  grid-distributed  flow  field  is  included  in  the  stack  analysis.  Combined  with  the  coordinate  change  preprocessing  and  analog 
technique,  neural  network  is  used  to  treat  the  MAP  of  compressor  and  turbine  in  the  air  system.  Three  kinds  of  air  system  topologies,  the  pure 
screw  compressor,  serial  booster  and  exhaust  expander  are  analyzed  in  this  article.  A  real-code  genetic  algorithm  is  programmed  to  obtain 
the  global  optimum  air  stoichiometric  ratio  and  the  cathode  outlet  pressure.  It  is  shown  that  the  serial  booster  and  expander  with  the  help  of 
exhaust  recycling,  can  improve  more  than  3%  in  the  FCS  efficiency  comparing  to  the  pure  screw  compressor.  As  the  net  power  increases,  the 
optimum  cathode  outlet  pressure  keeps  rising  and  the  air  stoichiometry  takes  on  the  concave  trajectory.  The  working  zone  of  the  proportional 
valve  is  also  discussed.  This  presented  work  is  helpful  to  the  design  of  the  air  system  in  fuel  cell  system.  The  steady-state  optimum  can  also 
be  used  in  the  dynamic  control. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Fuel  cells  provide  an  environmentally  friendly  high  effi¬ 
ciency  power  source  that  is  not  limited  by  the  Carot  effi¬ 
ciency.  The  polymer  exchange  membrane  fuel  cell  (PEMFC) 
is  considered  to  be  the  most  promising  candidate  for  electric 
vehicles  by  virtue  of  its  high  power  density,  zero  pollution, 
low  operating  temperature,  quick  start-up  capability  and  long 
lifetime.  PEMFC  can  also  be  used  in  distributed  power  sys¬ 
tems,  submarines,  and  aerospace  applications. 

The  properties  of  the  stacks  and  the  air  supply  subsystem 
are  essential  to  the  fuel  cell  system  performance.  The  system 
topology  and  operating  conditions  are  the  dominant  factors 
should  be  paid  attention  to.  Increasing  the  operating  pressure 
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and  air  flow  rate  improves  the  stack  power,  but  aggravates  the 
auxiliary  consumption  on  the  other  hand.  So  there  must  be 
a  balancing  point  in  the  operating  condition.  It  is  necessary 
to  analyze  the  trade-off  between  stacks  and  the  air  supply 
system  by  the  method  of  modeling  and  optimization. 

In  general,  the  stack  models  are  divided  into  the  empirical 
and  analytical  ones.  Kim  et  al.  [1]  introduced  the  semi- 
empirical  equations  to  describe  the  activation,  ohm  and  con¬ 
centration  overpotential.  Based  on  Amphlett’s  SSEM  [2], 
Mann’s  GSSEM  [3]  can  get  reasonable  fuel  cell  polarization 
by  considering  the  water  content  in  the  membrane.  Although 
the  empirical  and  semi-empirical  models  can  interpret  the 
polarization  and  the  effects  of  the  pressure,  they  both  almost 
failed  to  explain  the  influence  of  the  air  stoichiometry.  The 
water  transportation  is  also  omitted,  which  is  important  to 
the  water  and  thermal  management.  The  work  of  Bernardi 
and  Verbrugge  [4]  and  Springer  et  al.  [5]  are  the  milestone 


C.  Bao  et  al.  /  Journal  of  Power  Sources  156  (2006)  232-243 


233 


Nomenclature 

a  vapor  activity 

A  active  cell  area  (cm2) 

cref  reference  oxygen  concentration  (mol  m-3) 

Do  mutual  diffusion  coefficient  (m2  s-1) 

EW  equivalent  weight  of  polymer  ( 1 . 1  Kg  m-3 ) 

F  Faraday  constant  (96487  C  mol- 1 ) 

gear  transmission  ratio  in  Mode  III 

I  current  density  (A  cm-2) 

^drop  pressure  drop  coefficient  (bar  s  Kg- 1 ) 

m  mass  flow  rate  (Kg  s-1) 

M  molecular  weight  (Kg  mol  - 1 ) 

N  revolution  (rpm)  or  flux  (mol  s  - 1 ) 

ACe]i  number  of  cells  in  stacks 

N  compressor  speed  factor  (rpm  K-1/2) 
op_  optimum  value 

p  pressure  (bar) 

P  power  (kW) 

R  universal  gas  constant  (8.3 143  J  mol-1  K) 

SR  stoichiometric  ratio 

T  temperature  (K) 

wo  mass  fraction  of  water  within  the  membrane 

v  species  molar  fraction  or  neural  network  inputs 

y  mass  fraction  or  target  values  in  neural  network 

training 

z  length  in  GDL 

Greek  letters 

a  net  water  transport  coefficient 

s  porosity  of  GDL 

0  efficiency  or  overpotential  (V) 

k  electrical  conductivity  of  membrane 

(£2-1  m  1 ) 

X  water  content  in  membrane 

(mol  H2 O/equivalent  SO3-1) 
fi  dynamic  viscosity  (Kgm-1  s-1)  or  chemical 

potential  (J  mol-1) 

§  electro-osmotic  drag  coefficient  of  water  in 

membrane 

jt  pressure  ratio 

p  density  (Kgm-3) 

0  mass  flow  rate  factor  of  compressor 

(Kg  s-1  K1/2  bar-1)  or  potential  of 
electrolyte  (V) 

Superscripts  and  subscripts 
a  anode 

aux  auxiliary 

c  cathode  or  compressor 

H2  hydrogen 

in  inlet 

m  mechanical 

out  outlet 

O2  oxygen 


surge  surge  point  for  centrifugal  compressor 

t  turbine 

w  water  vapor 


for  the  analytical  PEM  fuel  cell  models.  Based  on  the  con¬ 
centrated  solution  theory,  Fuller  and  Newman  [6]  built  a 
more  compact  model  to  describe  the  three  water  transporta¬ 
tion  mechanisms.  Cheng  et  al.  [7]  analyzed  the  influence  of 
structural  parameters  and  operating  conditions  on  properties 
of  the  fuel  cell.  With  the  representative  work  of  Natarajan 
and  Nguyen  [8]  and  Wang  et  al.  [9],  the  analytical  fuel  cell 
models  have  been  developed  into  three-dimensional,  two- 
phase,  non-isothermal,  transient  ones  nowadays.  However, 
one-dimensional  models  are  more  suitable  for  the  system- 
level  optimization. 

Many  papers  have  been  written  on  pressurization  effects. 
Wiartalla  et  al.  [10]  compared  many  kinds  of  compres¬ 
sor/expander  technologies  and  their  suitability  for  use  in  fuel 
cell  systems.  Cunningham  et  al.  [11]  analyzed  the  differ¬ 
ence  between  high-pressure  and  low-pressure  systems  and 
introduced  a  flexible  air  supply  model  [12].  Friedman  et  al. 

[13]  analyzed  the  overall  balancing  demands  in  an  indirect 
methanol  PEM  fuel  cell  system.  And  Thirumalai  and  White 

[14]  simulated  the  steady-state  operation  of  compressor  with 
the  constant  oxygen  stoichiometry  or  the  constant  compres¬ 
sor  speed.  In  summary,  the  topology  and  operating  conditions 
are  the  two  major  topics  in  performing  steady-state  simula¬ 
tion  on  the  air  supply  system,  which  are  just  the  focus  of  this 
paper. 

An  analytical  model  is  adopted  in  this  article  to  simulate 
the  fuel  cell.  Taking  the  number  of  cells  into  consideration 
and  incorporating  the  simple  channel-groove  pressure  drop 
model,  the  single  cell  model  is  extended  into  a  stack  model. 
The  neural  network  technology  is  used  to  interpolate  and 
extrapolate  the  MAP  of  compressor  and  turbine.  Three  kinds 
of  topologies  of  air  supply  systems  are  discussed.  And  the 
real-code  genetic  algorithm  is  used  to  optimize  the  air  stoi¬ 
chiometric  ratio  and  the  cathode  outlet  pressure. 

2.  Systematic  topology 

System  topology  is  a  dominant  factor  in  determination 
of  the  operating  mode  of  a  fuel  cell  system.  Fig.  1  shows 
three  kinds  of  air  supply  subsystems  in  a  high-pressure  type 
application.  Mode  I  is  a  single-stage  scheme  with  a  sole 
screw  compressor  shown  in  Fig.  1(a).  Because  the  screw  com¬ 
pressor  has  a  wide  operating  zone  of  high  efficiency  and  is 
insensitive  to  the  presence  of  liquid  water,  this  kind  of  air 
system  is  widely  used  in  the  current  applications.  Modes  II 
and  III  both  utilize  the  energy  of  exhaust,  which  is  similar  to 
the  manner  in  the  internal  combustion  engine  (ICE).  Mode  II 
is  a  serial  layout  as  shown  in  Fig.  1(b).  Exhaust  gas  propels 
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Fig.  1.  Schematic  of  three  topologies  of  air  supply  subsystems,  where  PV, 
T,  C  represents  the  proportional  valve,  turbine  and  compressor,  respectively. 
The  arrowheaded  solid  line  is  for  the  intake  air  and  the  arrowheaded  dashed 
line  is  for  the  exhaust  gas:  (a)  is  Mode  I  of  single-level  screw  compressor, 
(b)  is  Mode  II  of  serial  booster  and  (c)  is  Mode  III  of  expander. 

the  centrifugal  compressor/turbine  to  supply  the  first-stage 
pressure  ratio  (71^1),  while  the  screw  compressor  obtains  the 
second-stage  pressure  ratio  (71^2).  Mode  III  is  the  so-called 
expander  as  depicted  in  Fig.  1(c).  An  auxiliary  motor  drives 
the  centrifugal  compressor  to  obtain  the  required  flow  rate 
and  pressure,  and  the  turbine  recycles  partial  exhaust  energy. 
The  proportional  valve  (PV)  is  optional  in  Modes  II  and  III. 


3.  Fuel  cell  system  modeling 

3.1.  Fuel  cell  model 


the  hydrogen  and  air  stoichiometric  ratio  and  the  net  water 
transport  coefficient  [5]: 


r°ut  _ 
Aw,a 


SR;l.Ca  -  2a(l  -  *“  „) 


w,a- 


in 

w,a 


-2a(l-  x™  a)  +  SRa  -  1 


T°ut  _ 
Aw,c 


SRC<C  +  0.42(1  +  2a) (1  -  <c) 


W,C' 


SRC+  0.21(1  +4a)(l-je“  ) 


xOU,  =  0.21(SRc-1)(1-<c) 

°2  SRC  +  0.21(1  +  4a)(l  -  c) 

In  order  to  reduce  the  two-dimensional  or  pseudo  two- 
dimensional  model  to  a  one-dimensional  model,  the  loga¬ 
rithmic  average  of  oxygen  concentration  and  the  arithmetic 
average  of  other  species  molar  fraction  between  the  inlet  and 
outlet  are  adopted  [2].  The  averaged  values  are  used  as  the 
boundary  condition  of  gas  diffusion  layer.  The  logarithmic 
average  of  oxygen  concentration  is  similar  to  the  logarith¬ 
mic  average  of  temperature  difference  in  the  heat  exchanger. 
This  treatment  is  the  key  to  reflect  the  effects  of  the  air  stoi¬ 
chiometry  as  shown  in  Fig.  2(a),  which  compares  the  voltage 
variation  with  the  increasing  air  stoichiometric  ratio  (SR)  pre¬ 
dicted  by  several  empirical  models  and  the  analytical  model 
in  this  paper.  Natarajan  and  Nguyen  [8]  proved  that  using 
the  inlet  molar  fractions  leads  to  notable  discrepancy.  As 
Fig.  2(a)  shows,  most  of  the  empirical  equations  underes¬ 
timate  the  influence  of  the  SR  on  the  fuel  cell  performance. 
When  these  empirical  models  are  used  in  the  optimization  of 
SR,  the  optimum  air  stoichiometric  ratio  is  almost  equal  to 
unity  because  the  small  gain  in  the  stack  power  cannot  offset 
the  increment  of  auxiliary  consumption. 

The  multicomponent  mass  transfer  in  the  gas  diffusion 
layer  (GDL)  is  described  by  Stefan-Maxwell  equation.  The 
effective  diffusion  layer  thickness  in  one-dimensional  model 
is  approximated  to  the  original  size  multiplied  by  1/0.6,  thus 
the  effects  of  serpentine  ribs  can  be  neglected  [15]: 


The  fuel  cell  model  in  this  article  is  a  modified  model 
based  on  antecessors’  work.  The  complexity  of  the  fuel  cell 
was  simplified  with  the  following  assumptions: 

1.  The  gas  diffusion  layers  and  membrane  are  all  isotropic. 
And  the  catalyst  layers  are  thin  infinitely. 

2.  Only  steady  state  is  considered. 

3.  The  cell  temperature  and  the  total  pressure  remain  con¬ 
stant. 

4.  All  species  are  in  the  gaseous  phase  and  the  phase  change 
of  water  was  not  included.  When  the  vapor  pressure  is  over 
than  the  saturated  pressure,  the  dilution  of  oxygen  can  be 
approximated  to  be  a  liquid’s  obstacle  to  mass  transfer  to 
some  extent. 

5.  Only  the  direction  normal  to  the  catalyst  layer  is  consid¬ 
ered  here. 

The  molar  fraction  of  each  species  at  the  outlet  of  flow 
channel  can  be  calculated  from  the  mass  conservation  given 


The  species  fluxes  are  related  to  current  density  and  net 
water  transport  coefficient: 


I  I 

^w,a  =  Oi  —  ,  Nh2  =  — , 


-(1  +  2a) 


I 

IF 


The  concentration  solution  theory  is  used  in  the  mem¬ 
brane  modeling  [6]  due  to  better  numeric  characteristics  than 
Springer’s  model  [5]  and  consideration  on  the  effects  of  pres¬ 
sure  difference  naturally: 


I  I 

a—  =  § 


F 


F 


pDo  dwo 

- — - -  VA 

(1  —  ujq)  Mh2o  dA 
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§  d/XH20 

F  dX 


VA 


1.98  +  0.0324A\ 
1  +  0.0648A  ) 


where  p,  wo.  Do ,  /c,  §  and  /xh2o  all  depend  on  the  water 
content  of  the  membrane  X: 


k  =  (0.5139A 


0.326)  exp 


AMH2o 


AMh2o  +  EW 

RTWa 
V/xh2o  = - 

a 


(11) 

(12) 

(13) 


Dq  =  3.5  x  10  6  exp 


-2436 

T 


X 

14 


1 


1 


|_(0.35A)4  +  (1.47)4  J 


-1/4 


(9) 

(10) 


The  water  content  in  both  sides  of  the  membrane  is  related 
to  vapor  activity  a: 


'A  =  0.043  +  17.81 a 
—39.85 a2  +  36.0  a3 
X  =  14  +  2.8g-1(a  —  1) 
„  X  =  16.8 


fora  <  1 

for  1  <  a  <  1  +  g 
for  a  >  1  +  g 


(14) 


Here,  the  value  of  g  is  30  as  Hsing  [  1 6]  adopted  instead  of  2 
used  in  the  Springer’s  model  [5].  Assuming  the  value  of  a ,  we 
can  obtain  two  values  of  the  water  content  in  the  cathodic  end 
of  the  membrane.  One  is  from  equations  of  the  anode  GDL 
and  membrane  and  the  other  is  from  equations  for  cathode 
GDL.  With  the  boundary  condition  that  both  water  contents 
should  be  identical,  the  Bvp4c  algorithm  in  MATLAB  shows 
better  numeric  characteristic  than  Newton-Raphson  algo¬ 
rithm  to  solve  for  a.  The  open  circuit  voltage  is  from  the 
Nernst  equation,  and  the  ohmic  overpotential  can  be  obtained 
by  integrating  Eq.  (7).  The  simple  Tafel  expression  is  used  to 
calculate  the  cathode  activation  loss: 


RT  f  I  p  x o2\ 
VaCt  ~  0.5 F  n  (/o  RT  ever ) 


(15) 


Fig.  2.  (a)  Comparison  of  the  influence  of  the  air  stoichiometry  on  the  fuel 
cell  voltage  increment  between  several  empirical  models  and  the  analytical 
model  in  this  paper  as  the  current  density  is  0.5  A  cm-2,  (b)  Comparison 
of  the  predicted  polarization  curve  with  experimental  results  as  the  air  sto¬ 
ichiometric  ratio  is  2.5.  The  cell  temperature  is  343  K,  hydrogen  utilization 
is  80%,  fuel  and  air  are  both  saturated,  the  total  pressure  is  3  bar  and  there 
is  no  pressure  difference  between  the  anode  and  the  cathode.  The  effective 
thickness  of  GDL  is  0.056  cm,  porosity  of  GDL  is  0.35,  and  the  open  circuit 
voltage  is  0.975  V  which  is  from  experiment  instead  of  Nernst  equation. 


The  predicted  polarization  curve  is  compared  with  exper¬ 
imental  results  in  Fig.  2(b)  with  reasonable  agreements  at 
low  and  moderate  electrical  loads.  This  model  can  also  be 
used  in  water  and  thermal  management  of  fuel  cell  system 
[17]. 

3.2.  Pressure  drop  model 

In  addition  to  the  simple  extension  from  a  single  cell 
to  stacks,  the  pressure  drop  model  of  air  flow  field  is  also 
included.  Fig.  3  shows  the  grid-distributed  prototype  of 
the  pressure  drop  model.  For  simplicity  the  ribbed  car¬ 
bon  plates  is  regarded  to  form  a  rectilinear  network  of 
grooves  [18]  and  channels  after  making  some  assumptions  as 
follows: 

1.  The  pressure  drop  in  the  air  intake  manifold  is  not 
considered  and  the  flow  rate  is  uniform  in  each  single 
cell. 

2.  The  fluid  is  under  a  fully-developed  laminar  flow,  of  which 
the  expansion  and  contraction  losses  of  fluid  occurring  due 
to  the  change  in  flow  directions  are  neglected. 

3.  The  current-density  distribution  along  the  length  and 
breadth  of  a  single  fuel  cell  is  assumed  to  be 
uniform. 
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Fig.  3.  The  prototype  of  the  grid-distributed  pressure  drop  model  in  the  cathode.  The  design  parameters  are  the  width  of  the  channel  ( H)  and  the  grooves  (V) 
and  depth  of  the  groove  (DP),  distance  between  channels  (DH)  and  grooves  (DV),  the  electrode  and  flow-field  width  (W)  and  length  (L),  and  the  width  of  inlet 
(VFi)  and  outlet  (Wo).  And  the  sizes  are:  L  =  40  cm,  W=  10  cm,  H- 0.25  cm,  V  =  0.025  cm,  DH  =  1  cm,  DV  =  0.4  cm,  DP  =  0.07  cm,  W\  =  1  cm  and  WQ  =  1  cm. 


For  laminar  flow  in  a  rectangular  duct,  the  relationship 
between  pressure-drop  and  flow-rate  is 


q  = 


4  ba3  —  A  p 
3  /J  A  l 


oo 


1  192a  tanh(inb/2a) 


7 r3b 


i=  1,3,5 


V 


q 


=  —kAp 


(16) 


where  a  and  b  are  the  half  width  and  half  depth  of  the  duct,  q 
the  volume  flow  rate  of  the  reactant  gas  and  Al  is  the  length 

of  the  duct. 

For  a  control  volume  at  the  node  P(iJ)  as  shown  in  Fig.  3, 
the  mass  conservation  equation  is  [19]: 

(/w  -  /e)  +  Un  ~  Js)  =  Sij  (17) 

where  Sq  is  the  consumption  of  the  reactant  in  the  control 
volume  and  /w  the  mass  flow  rate  related  to  the  node  W(i—  1  j) 
in  the  west: 

7 w  —  PwhiPw  P?)H  (18) 


The  density  in  the  boundary  of  control  volume  is  the  harmonic 
average  of  the  densities  at  two  adjacent  nodes  [19]: 

J_  =  (Ax)w+/Ax  +  (Ax)w_/Ax 
Pw  Pw  PE 

Eqs.  (16)— ( 19),  with  the  boundary  conditions  determined  by 
the  flow  through  the  cell,  form  a  set  of  algebraic  equations 
and  are  solved  using  TDMA  algorithm.  In  this  design  scale 
as  shown  in  Fig.  3,  if  the  current  density  is  0.4  A  cm-2,  air 
stoichiometric  ratio  is  2  and  the  inlet  pressure  is  2  bar,  and 
if  the  air  is  the  saturated  in  the  stack  operating  tempera¬ 
ture,  the  pressure  drop  coefficient  in  the  cathode  is  almost 


0.94  bar  s  Kg-1  without  considering  the  block  of  condensed 
water.  The  arithmetic  average  of  the  inlet  and  outlet  pressure 
is  used  as  the  total  pressure  in  the  single  cell  simulation. 

3.3.  Compressor/turbine  model 

In  general,  the  modeling  of  the  compressor  and  turbine  is 
based  on  their  MAP.  Here  neural  network  technology  is  used 
to  interpolate  and  extrapolate  the  MAP  because  of  its  good 
generalization  property. 

3.3.1.  Neural  network  architecture 

Artificial  neural  networks  are  models  of  information  pro¬ 
cessing  inspired  by  the  biology  of  the  brain.  The  basic  unit 
is  the  artificial  neuron  which  receives  numerical  informa¬ 
tion  from  input  nodes,  processes  it  internally,  and  gives  out 
outputs.  Networks  with  biases,  a  sigmoid  layer,  and  a  lin¬ 
ear  output  layer  are  capable  of  approximating  any  function 
with  a  finite  number  of  discontinuities.  Fig.  4  shows  the  neu¬ 
ral  network  architecture  in  the  compressors’  model  with  a 
hidden  layer  and  an  output  layer.  It  is  a  typical  topology  of 
feed-forward  multilayer  perceptron. 

At  first,  we  initialize  the  layer’s  weights  and  biases  accord¬ 
ing  to  the  Nguyen-Widrow  algorithm  which  chooses  values 
in  order  to  distribute  the  active  region  of  each  neuron  in  the 
layer  evenly  across  the  layer’s  input  space.  The  inputs  to  any 
processing  element  are  multiplied  by  connection  weights  and 
summed  in  sequence: 

neti  =  IW  x  +  bi, 


net2  =  LW  Oi  +  b2 


(20) 
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Fig.  4.  The  neural  network  architecture  of  compressor/turbine  model  with  a  hidden  layer  and  an  output  layer.  IW  is  the  weight  matrix  between  the  input  layer 
and  the  hidden  layer,  LW  is  the  weight  matrix  between  the  hidden  layer  and  the  output  layer,  bi  and  b2  are  the  biases  matrix  of  the  hidden  layer  and  the 
output  layer,  S\  the  number  of  neurons  in  the  hidden  layer.  The  transfer  functions  in  the  hidden  layer  and  output  layer  are  hyperbolic  tangent  and  linear  type, 
respectively. 


and  then  the  transfer  functions  are  used  to  calculate  a  layer’s 
output  from  its  net  input: 


Oi  =  /i(neti)  = 


2 

1  +  exp(— 2  •  neti) 


y  =  02  =  /2(net2)  =  net2 


(21) 


When  all  the  samples  have  been  calculated,  the  error  can  be 
obtained  by  comparing  the  target  values  with  the  network 
outputs: 


(22) 


Training  the  network  is  to  minimize  the  above  quadratic 
cost  function  by  adjusting  the  weights  and  biases.  Back  prop¬ 
agation  algorithm  is  the  most  commonly  used  optimization 
method  in  training  neural  networks.  Levenberg-Marquardt 
back  propagation  algorithm  is  used  in  this  paper  which  is  very 
suitable  to  the  optimization  of  quadratic  questions  because 
of  its  second-order  accuracy  and  no  need  to  calculating  the 
Hessian  matrix.  Given  two  of  the  three  variables,  compres¬ 
sor/turbine  speed  factor  (N  =  N/^/T),  pressure  ratio  (jtk)  or 
mass  flow  rate  factor  (0  =  m\ff  / p)  as  the  network  inputs, 
the  other  variable  can  be  obtained  from  the  trained  network. 
And  the  relationship  between  the  compressor  efficiency  (r/) 
and  relevant  parameters  0  and  it rj  =/(0,7T£),  can  be  obtained 
using  another  similar  network. 


3.3.2.  Preprocessing 

In  order  to  have  the  network  trained  well,  some  prepro¬ 
cessing,  such  as  normalization  is  inevitable.  Moraal  [20] 
compared  four  methods  treating  the  MAP  of  the  centrifu¬ 
gal  compressor/turbine,  and  found  the  neural  network  has 
the  best  performance.  But  he  failed  in  training  the  network  to 
approximate  the  nonlinear  relationship  0  =  /( 7z>,  A)  instead 
of  n k  =  /(0,  N).  That  is  because  the  compressor  speed  con¬ 
tours  shown  in  Fig.  5(a)  are  almost  parallel  to  the  axis  of 
flow  rate  (that  is  normal  to  the  axis  of  pressure  ratio)  in  the 
area  close  to  the  surge  line.  So  a  small  disturbance  of  the 
pressure  ratio  will  lead  to  a  large  variation  of  the  flow  rate, 


(b)  Mass  flow  factor  /Kg.s_1.K1/2.bar-1 

Fig.  5.  The  coordinate  transform  in  the  preprocessing  of  centrifugal  com¬ 
pressor  MAP  trained  by  neural  network:  (a)  shows  the  original  speed 
contours  and  (b)  shows  the  speed  contours  with  the  clockwise  rotated  coor¬ 
dinates. 
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which  causes  the  network  collapse.  In  this  paper,  coordinate 
transform  is  used  to  circumvent  this  problem. 


COS0 

sin# 

COS  6 

— sin# 


—sin  6 

0 

COS  6 

sin# 

~<Pr 

COS  6 

K 

0 


Jtk 


(23) 


where  0  is  the  coordinate  transform  angle  and  its  value  is  cho¬ 
sen  as  45°.  Through  this  treatment,  the  MAP  of  compressor 
looks  like  that  shown  in  Fig.  5(b),  and  the  neural  network  can 
have  wonderful  performance  in  the  case  of  0  =  /(7T&,  N). 

Another  method  for  improving  generalization  is  called 
early  stopping.  In  this  technique  the  available  data  are  divided 
into  three  subsets.  The  first  subset  is  the  training  set,  which 
is  used  for  computing  the  gradient  and  updating  the  network 
weights  and  biases.  The  second  subset  is  the  validation  set. 
The  error  on  the  validation  set  is  monitored  during  the  training 
process.  The  validation  error  will  normally  decrease  during 
the  initial  phase  of  training,  as  does  the  training  set  error. 
However,  when  the  network  begins  to  overfit  the  data,  the 
error  on  the  validation  set  will  typically  begin  to  rise.  When 
the  validation  error  increases  for  a  specified  number  of  itera¬ 
tions,  the  training  is  stopped,  and  the  weights  and  biases  at  the 
minimum  of  the  validation  error  are  returned.  Fig.  6  shows 
the  treatment  of  efficiency  MAP  of  the  screw,  centrifugal 
compressor  and  turbine,  rj  =/(</>, 77^)  with  good  interpolation 
and  extrapolation  accuracy  by  a  BP  network  with  7,  7  and  9 
neurons  in  the  hidden  layer,  respectively.  The  speed  contours 
can  be  easily  approximated  by  a  neural  network  because  of 
its  smooth  tendency. 


3.3.3.  Size  match 

In  this  paper,  all  the  MAP  of  screw,  centrifugal  compressor 
and  turbine  are  taken  from  reference  [21]  because  of  lack¬ 
ing  experimental  data.  The  original  charts  are  typical  for  a 
250  kW  stack,  so  the  sizes  of  stack  and  compressors  must  be 
matched  using  analog  technique  when  a  stack  with  different 
power  is  concerned.  According  to  the  analog  law,  two  com¬ 
pressors  are  the  analogue  if  they  share  common  Reynolds 
number  and  Mach  number.  Given  the  geometric  ratio  of  two 
analogical  compressors  as  r,  the  analog  characteristics  can  be 
obtained  as  follows: 

D 

—  =  r  =>  7tk  =  7tk,  t?  =  , 

N  1  P  o 

—  =  -,  —  =  r 

N'  r  P’ 


4.  Optimization  and  discussions 

4.1.  Optimization  problem 


For  a  fuel  cell  system,  the  optimization  target  is  to  obtain 
the  maximum  system  efficiency.  The  air  stoichiometric  ratio 
(SR)  and  the  cathode  outlet  pressure  (p0 ut)  are  the  two  most 
important  variables  influencing  the  trade-off  between  the 
stack  and  the  air  supply  system.  With  a  constant  hydrogen 
stoichiometric  ratio,  the  equivalent  objective  is  to  obtain  the 
maximum  net  power.  The  optimization  problem  in  this  paper 
can  be  described  as  follows: 

maximize  /(SR,  p0 ut)  =  Pfcs  =  Pstack  P motor  Paux 

(25) 


where,  the  motor  auxiliary  power  Pmotor  is  the  driving  motor 
consumption  of  the  screw  compressor  in  Modes  I  and  II,  and 
the  auxiliary  motor  consumption  of  the  expander  in  Mode  III. 
The  power  Paux  used  by  all  other  auxiliary  devices  means  the 
power  sum  of  pumps  and  fans  and  so  on,  which  is  chosen  as 
1  kW  in  this  paper. 

For  Mode  I,  the  optimization  is  an  unconstraint  problem 
within  the  searching  range.  But  for  Mode  II,  there  are  con¬ 
straints  as  follows: 


< 


Pin  —  Pout  T  Kft 
Pc  —  P 01m 


rop 


m 


y  air 


1.01  3tC  fc\TC 


,  Pout  >  TTt  (PV  works)  or  pout  =  Tt{  (without  PV) 


(26) 


For  Mode  III,  it  is  subject  to  the  constraints: 

f  m 

Pin  —  Pout  T  ^drop  <  Psurge 

)7ur 
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Nc  =  M  gear 


,  out  >  7Tt  (PV  works)  or  pout  =  Ttt  (without  PV) 


Where, 

I A  SR  N cep  Mair 

m  =  - 

4F0.21 


(27) 


(28) 


Note  that  there  is  a  transmission  ratio  between  the  compressor 
and  the  turbine  in  Mode  III  in  order  to  widen  the  united  work¬ 
ing  range.  It  is  found  that  the  expander  can  just  start  to  work 
effectively  when  the  power  is  over  50  kW  without  a  transmis¬ 
sion  ratio.  One  of  the  reasons  behind  may  be  that  the  original 
combination  of  compressor/turbine  is  not  suitable  to  work  as 
an  expander  in  this  relatively  smaller  application.  Table  1  lists 
the  stack  and  compressors  sizes  in  the  three  modes.  The  rated 
net  power  of  fuel  cell  system  is  60  kW.  Here  the  analog  ratio 
of  compressors  represents  the  geometric  ratio  to  the  original 
size,  which  is  chosen  by  balancing  between  high  FCS  net 
power  and  wide  working  range  of  high  efficiency. 


Such  relationship  is  also  suitable  for  the  screw  compressor 
and  the  turbine.  It  is  obvious  that  there  should  be  different 
analog  ratios  of  compressors/turbine  in  different  systematic 
topologies  which  will  be  discussed  later. 


4.2.  Real-code  genetic  algorithm 

In  this  article,  a  real-code  genetic  algorithm  (GA)  is  used 
to  obtain  the  global  optimum  SR  and  poui  as  the  current  den- 
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Fig.  6.  The  efficiency  MAP  of:  (a)  screw  compressor,  (b)  centrifugal  compressor,  (c)  turbine  processed  by  a  BP  neural  network  with  7,  7  and  9  neurons  in  the 
hidden  layer,  respectively.  Dashed  lines  are  the  original  data,  solid  contours  are  the  trained  interpolated  results  and  dotted  contours  represent  the  extrapolated 
network  outputs. 


sity  varies.  GA  is  a  search  procedure  based  on  the  mechanism 
of  natural  selection.  It  uses  operations  found  in  natural  genet¬ 
ics  to  guide  itself  through  the  paths  in  the  search  space.  GA 
requires  the  problem  of  maximization  (or  minimization)  to 
be  stated  in  the  form  of  a  cost  function.  In  GA,  a  set  of  vari- 

Table  1 


Sizes  of  stacks  and  compressors  in  the  three  modes 


Active  cell  area  (cm2) 

400 

Number  of  cells  in  one  stack 

330 

Number  of  serial  stacks 

1 

Number  of  parallel  stacks 

2 

Membrane  type 

Nafion  112 

Analog  ratio  of  screw  in  Mode  I 

0.7071 

Analog  ratio  of  screw  in  Mode  II 

0.7071 

Analog  ratio  of  compressor  in  Mode  II 

0.7156 

Analog  ratio  of  turbine  in  Mode  II 

0.9090 

Analog  ratio  of  compressor  in  Mode  III 

0.6794 

Analog  ratio  of  turbine  in  Mode  III 

0.9258 

Transmission  ratio  in  Mode  III 

0.7 

ables  is  encoded  into  a  string,  analogous  to  a  chromosome 
in  nature.  To  determine  how  well  a  chromosome  solves  the 
problem,  it  is  first  broken  down  into  the  individual  substrings 
which  represent  each  variable  and  these  values  are  then  used 
to  evaluate  the  cost  function,  yielding  fitness.  GA  selects  par¬ 
ent  from  a  pool  of  population  according  to  the  basic  criteria 
of  “survival  of  the  fittest”.  The  repopulation  of  the  next  gen¬ 
eration  is  done  using  three  methods:  reproduction,  crossover, 
and  mutation. 

Compared  to  the  traditional  binary-code  GA,  a  real-code 
GA  has  clearer  physical  meaning  in  the  chromosome  encod¬ 
ing  and  comprises  some  steps  as  follows: 

•  Step  1 :  generate  a  cluster  with  30  random  chromosomes 
as  [SR,  pout]  and  define  the  cost  function  as  the  FCS  effi¬ 
ciency. 

•  Step  2:  copy  and  cross  chromosomes  by  twos.  The 
crossover  of  two  chromosomes  in  the  real-code  GA  is 


Cathode  Outlet  Pressure  (bar) 
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defined  by 

1 ^ik=®ik(\  fik)  +  Ojkfiki  djk  ~  &jk(  1  fik)  T  dikfik 

(29) 

where  is  the  kth  variable  in  the  ith  chromosome  of  a 
cluster.  To  interpolate  two  offspring  as  the  random  number 
Pk  is  within  [0,  1]  while  extrapolate  two  offspring  as  fik  is 
within  [—1,0]  from  the  parent. 

•  Step  3:  calculate  the  fitness  of  90  chromosomes  and  select 
30  better  chromosomes. 

•  Step  4:  mutate  in  a  given  probability.  And  the  mutation  of 
a  chromosome  is  defined  by 

aik  =  aik  -  ( aik  -  afm)m  forO  <  &  <  0.5 
a'ik  =  aik  +  («“aX  “  aik)m k  for  0.5  <  ffc  <  1 

where  [a™m,  a™m]  is  the  range  of  Mi  variable,  and  rjk  and 
are  two  random  numbers. 


•  Step  5:  calculate  the  fitness  of  new  cluster  and  compete 
chromosomes  by  the  roulette  selection  which  ensures  the 
better  chromosome  has  big  probability  to  win. 

•  Step  6:  judge  whether  the  difference  between  maximum 
fitness  of  two  generations  satisfies  the  tolerance  and  there 
is  no  better  chromosome  within  the  five  continuous  gen¬ 
erations.  If  it  satisfies,  algorithm  converges  and  the  chro¬ 
mosome  with  maximum  fitness  is  the  solution,  else  go  to 
the  step  2. 

The  real-code  GA  is  validated  better  than  a  simplex  algo¬ 
rithm  to  get  global  optimum  in  this  study  case. 

4.3.  Results  and  discussions 

4.3.1.  Optimization 

Fig.  7  shows  the  searching  progresses  for  the  three  modes 
of  the  topology  of  the  air  system  with  the  parameters  listed  in 
Table  2  in  the  rated  condition  under  which  the  current  density 


Fig.  7.  Efficiency  Contours  of  the  fuel  cell  system  corresponding  to  the  air  stoichiometric  ratio  (SR)  and  the  cathode  outlet  pressure  (pou t)  in  the  rated  condition, 
which  means  the  current  density  is  0.4  A  cm-2  with  the  basic  operating  conditions  listed  in  Table  2;  (a)-(c)  is  for  Mode  I— III,  respectively. 
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Table  2 

Basic  operating  conditions 

Operating  cell  temperature  (K)  353 

Hydrogen  stoichiometry  1 .25 

Anode  inlet  relative  humidity  (%)  100 

Cathode  inlet  relative  humidity  (%)  100 

Pressure  difference  between  anode  and  cathode  (bar)  0 

Motor  efficiency  (%)  80 

Other  total  auxiliary  consumption  (kW)  1 


is  0.4  A  cm-2.  The  contours  of  the  optimizing  variables  show 
the  trade-off  between  stack  and  air  system  clearly.  With  the 
help  of  recycling  exhaust,  Modes  II  and  III  give  higher  effi¬ 
ciency  in  bigger  air  flow  rate,  which  also  benefit  to  water  and 
thermal  management.  The  lower  boundary  of  the  pressure 
obtained  for  Modes  I  and  II,  each  of  which  forms  a  curve¬ 
line  with  changeable  air  stoichiometric  ratio,  represents  flow 
resistance  in  the  turbine,  while  the  upper  boundary  of  the 
pressure  for  Mode  III  is  due  to  the  surge  limit  of  the  centrifu¬ 
gal  compressor. 

Fig.  8  shows  the  optimum  air  stoichiometry  (op_SR)  and 
cathode  outlet  pressure  (op_pGut)  of  the  three  modes  of  air 
system  topology  as  current  density  changes  from  0.01  to 
0.7  A  cm-2  while  other  operating  conditions  are  the  same  as 
that  in  Table  2.  With  an  increasing  net  FCS  power,  op_pout 
keeps  rising  for  all  three  modes,  and  op_SR  reveals  the 
concave  trajectory  except  under  small  load  condition.  To 
avoid  the  numeric  problem  in  relation  to  the  flat  property  of 
screw  efficiency  contour,  the  searching  lower  limit  of  pout 
is  chosen  as  1.3  bar  which  is  reasonable  in  high-pressure 
system.  That  is  why  op_SR  keeps  ascending  in  the  case  of 
small  electric  loads  in  Modes  I  and  II.  If  there  is  no  such 
limit,  op_SR  in  Modes  I  and  II  will  keep  descending  and 
will  be  not  less  than  5  when  the  current  density  is  very 
small,  while  op_pGut  will  be  almost  equal  to  the  atmosphere 
pressure  in  small  load  zone.  Because  the  compressor 
surge  limit  masks  the  limit  of  1.3  bar  naturally,  the  air 
stoichiometry  descends  in  Mode  III  when  the  net  power  is 
low. 

As  shown  in  Fig.  8(a),  the  optimum  air  stoichiometric 
ratios  in  the  Modes  II  and  III  are  much  higher  than  that 
in  Mode  I  because  of  the  exhaust  recycling.  The  corre¬ 
sponding  net  power  at  the  turning  point  of  op_SR  of  Mode 
I  is  the  highest  while  that  of  Mode  III  the  lowest.  The 
turning  points  of  op_SR  in  Modes  II  and  III  are  the  united 
working  timing  of  compressor/turbine,  so  the  expander 
starts  to  work  effectively  earlier  than  the  serial  booster 
does. 

Fig.  9(a)-(f)  shows  the  optimization  trajectory  of  the 
screw  compressor,  centrifugal  compressor  and  turbine  in  all 
three  modes.  In  Mode  I,  the  optimum  line  is  in  the  normal 
direction  of  the  screw  efficiency  contours.  In  Mode  II,  the 
optimum  condition  dominated  by  turbine,  is  along  the  normal 
direction  of  turbine  efficiency  contours.  Because  of  the  low 
energy  of  the  exhaust  gas,  the  first-stage  pressure  ratio  is 


Fig.  8.  The  optimum  air  stoichiomtric  ratio  (SR)  and  cathode  outlet  pres¬ 
sure  (/?out)  three  modes  related  to  the  FCS  net  power:  (a)  is  for  the  SR 
optimization  and  (b)  is  for  the  pressure  optimization. 

so  low  that  the  compressor  is  far  from  the  high  efficiency 
zone.  So  the  centrifugal  compressor  is  unsuitable  to  be  used 
in  the  serial  boost  mode  unless  the  slope  of  the  surge  line  is 
decreased.  In  the  expander  scheme,  the  optimum  trajectories 
as  shown  in  (Fig.  9(e)  and  (f))  are  balanced  between 
the  compressor  and  the  turbine.  When  the  compressor  is 
dominant,  the  optimum  solution  is  close  to  compressor  surge 
line,  or  it  is  along  the  normal  direction  of  turbine  efficiency 
contours. 

Fig.  9(g)  compares  the  fuel  cell  system  efficiencies 
in  the  three  modes.  With  the  help  of  exhaust  recycling, 
the  system  efficiency  can  be  elevated  more  than  3% 
relative  to  the  sole  screw  scheme.  Providing  the  prop¬ 
erty  of  the  centrifugal  compressor  changeable,  the  serial 
booster  maybe  given  higher  system  efficiency  than  the 
expander  does  because  the  compressor  efficiency  is  related 
to  the  pressure  ratio  exponentially  and  to  the  flow  rate 
linearly. 
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(a)  Mass  flow  factor  / Kg.s'kK^.bar'1 
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(b)  Mass  flow  factor  /Kg.s'1.K1/2.bar'1 


(c)  Mass  flow  factor  /Kg.s'kK^.bar'1  (d)  Mass  flow  factor  /Kg. s'1. K1/2. bar'1 


(e)  Mass  flow  factor  /Kg.s'VK^.bar'1  (f)  Mass  flow  factor  /Kg. s'1. K1'2. bar"1 


Fig.  9.  (a)-(f)  is  the  optimum  trajectories  of  the  compressors  and  turbine  in  three  modes:  (a)  is  for  the  screw  in  Mode  I;  (b)-(d)  is  for  the  centrifugal  compressor, 
turbine  and  screw  in  Mode  II,  respectively;  (e)-(f)  is  for  the  centrifugal  compressor  and  the  turbine  in  Mode  III.  The  dashed  line  in  (d)  and  (e)  means  that  the 
turbine  is  bypassed,  (g)  Shows  the  FCS  efficiencies  in  the  three  modes. 


4.3.2.  The  proportional  valve  working  zone 

In  Modes  II  and  III,  the  proportional  valve  may  not  work 
in  some  conditions.  Fig.  10  compares  the  optimum pQ ut  of  the 
expander  schemes  with  and  without  the  proportional  valve. 
Under  both  conditions,  op_p0ut  are  identical  in  the  case  of 
middle  electric  loads,  which  mean  the  proportional  valve  is 
fully  open  and  the  outlet  pressure  is  completely  controlled  by 


the  flow  rate.  When  the  electric  load  is  small,  the  turbine  does 
not  work  and  the  proportional  valve  must  be  there.  In  the  case 
of  high  electric  loads,  the  optimum  pressure  requires  that  the 
proportional  valve  keeps  partly  open.  If  the  variable  geome¬ 
try  turbine  (VGT)  is  adopted,  the  proportional  valve  can  be 
removed  and  the  characteristics  of  the  air  supply  system  will 
be  more  flexible. 
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Fig.  10.  The  optimum  cathode  outlet  pressure  in  the  expander  mode  with 
and  without  the  proportional  valve  as  the  current  density  varies  from  0.01 
to  0.7  A  cm~2,  while  other  operating  conditions  are  the  same  as  listed  in  the 
Table  2. 

5.  Conclusions 

The  contributions  in  this  paper  can  be  concluded  as  fol¬ 
lows: 

1.  An  analytical  fuel  cell  model  was  used  and  logarithmic 
mean  of  oxygen  pressure  was  adopted  to  avoid  underes¬ 
timating  the  influence  of  SR.  The  simple  channel-groove 
pressure  drop  model  was  included  in  the  stack  analysis. 

2.  Neural  network  was  used  to  treat  the  MAP  of  compressors, 
which  was  combined  with  the  coordinate  transfer  prepro¬ 
cessing.  And  analog  technology  was  adopted  to  match 
parameters  of  the  air  system. 

3 .  U sing  a  real-code  genetic  algorithm,  the  air  stoichiometric 
ratio  and  the  cathode  outlet  pressure  were  optimized  and 
compared  for  three  topologies  of  the  air  supply  systems. 

4.  Using  exhaust  recycling,  FCS  efficiency  can  be  improved 
over  3%  than  that  in  the  sole  screw  scheme.  The  working 
zone  of  the  proportional  valve  was  also  discussed. 

The  methodology  of  configuration  and  optimization  pre¬ 
sented  in  this  paper  is  valuable  for  design  and  analysis  of  the 
air  system  in  fuel  cell  system.  And  the  steady-state  results 
can  also  be  used  in  the  dynamic  control. 
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